      program rkamb72b
      implicit real*8 (a-h,o-z)
      parameter(np=3,mp=np+1,nplay=101,ndec=12,
     &          ftol=1.0d-10)
c
      external famoeb
      character*6 name(np)
      character*24 fdate
      dimension p(mp,np),x(np),y(mp),z(nplay,np)
      dimension is(nplay,ndec)
      common /mats/is
      common /iplay/id
c
      data name/'GAMMA','BETA2','BETA1'/
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c      do 5 i=1,nplay
c         write(*,'(i2,2x,12(i1,x))')i,(is(i,n),n=1,ndec)
c    5 continue
c      stop
c
c********************************************************************** 
      write(*,'('' ***** rkamb72b ******'',20x,a24,
     &//"    Ambiquity & Het Experiment")')fdate() 
      write(*,'(//"These are the amoeba parameter estimates:")')
      write(*,'(/" ID     Gamma         BETA2       BETA1      ",
     &"LL"/)')
c
      cumll=0.d0
      do 100 id=1,nplay
c
         x(1) = dlog(1.d0/0.3d0 - 1.d0)
         x(2) = dlog(1.d0/0.9d0 - 1.d0)
         x(3) = dlog(1.d0/0.1d0 - 1.d0)

c initialize starting vertices of the simplex. 
         call init(x,p,y,famoeb)
c 
c call the simplex routine. 
         ndim=np 
         iter=3000
c         iter=1
         call amoeba(p,y,mp,np,ndim,ftol,famoeb,iter) 
c 
c output results and stop. 
c         write(*,'(a,i5)')'iterations: ',iter 
c compute average values of vertices. 
         call avg(p,x,np,mp) 
         xllf=-famoeb(x) 
         cumll = cumll + xllf
         do 75 k=1,np
            z(id,k) = x(k)
   75    continue
c reverse parameter transformation
         x(1)=2.d1/(1.d0+dexp(x(1)))
         x(2)=1.d0/(1.d0+dexp(x(2)))
         x(3)=x(2)+1.d0/(1.d0+dexp(x(3)))
         if (x(3).gt.1.d0) x(3)=1.d0
c
         write(*,'(i3,2x,4g13.5)')id,(x(k),k=1,np),xllf
  100 continue
      write(*,'(/32x,g13.6)')cumll
c      call xpost(z)
      end
c***********************************************************************
      double precision function famoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=3,nplay=101,ndec=12)
      dimension x(np),px(ndec,2)
      dimension is(nplay,ndec)
      common /mats/is
      common /iplay/id
c
      gam=2.d1/(1.d0+dexp(x(1)))
      beta2=1.d0/(1.d0+dexp(x(2)))
      beta1=beta2 + 1.d0/(1.d0+dexp(x(3)))
      if (beta1 .gt. 1.d0) beta1=1.d0
c
      call fpx(gam,beta1,beta2,px)
c
      t=1.d0
      do 10 k=1,ndec
         j = is(id,k)
         t = t*px(k,j)
   10 continue
      famoeb = -dlog(t)
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta1,beta2,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      prj = (1.d0+beta2)/3.d0
c
      za = y1/2.d0
      zb = y1*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = y1/2.d0
      zb = y1*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = y1/2.d0
      zb = y2*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = y1/2.d0
      zb = y2*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = y1/2.d0
      zb = y3*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = y1/2.d0
      zb = y3*beta1/2.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta2*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta2*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta2*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*prj
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*prj
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*prj
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)
c
      return
      end
c***********************************************************************
      subroutine xpost(z)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,np=3)
      dimension px(ndec,2),z(nplay,np)
      dimension is(nplay,ndec)
      common /mats/is
      common /iplay/id
c
      do 100 id=1,nplay
         gam=2.d1/(1.d0+dexp(z(id,1)))
         beta2=1.d0/(1.d0+dexp(z(id,2)))
         beta1=beta2+1.d0/(1.d0+dexp(z(id,3)))
         if (beta1 .gt. 1.d0) beta1 = 1.d0
c
         call fpx(gam,beta1,beta2,px)

         t=1.d0
         do 10 k=1,ndec
            j = is(id,k)
            t = t*px(k,j)
   10    continue
         fl = dlog(t)
c
c         write(*,'(/"Subject ",i2,5x,g13.5)')id,fl
c         do 20 k=1,ndec
c            write(*,'(3x,i2,2g12.4)')k,(px(k,j),j=1,2)
c   20    continue
  100 continue
c
      return
      end
c***********************************************************************
      subroutine init(xint,p,y,famoeb)
      implicit real*8 (a-h,o-z) 
      parameter (np=3,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      do 2 k=1,np
c         if (xx(k).le.1.d-10) xx(k)=1.d-10
c         if (xx(k).ge.0.999) xx(k)=0.999
c         xxx(k)=dlog(1.d0/xx(k) - 1.d0)
c    2 continue
c
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
      y0=-famoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
      delta=0.5d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end

